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An event-based integration scheme for an integrate-and-fire neuron model with exponentially de- 
caying excitatory synaptic currents and double exponential inhibitory synaptic currents has recently 
been introduced by Carnevale and Hines. This integration scheme imposes non-physiological con- 
straints on the time constants of the synaptic currents it attempts to model which hamper the 
general applicability. This paper addresses this problem in two ways. First, we provide physical 
arguments to show why these constraints on the time constants can be relaxed. Second, we give a 
formal proof showing which constraints can be abolished. This proof rests on a generalization of the 
Carnevale-Hines lemma, which is a new tool for comparing double exponentials as they naturally 
occur in many cascaded decay systems including receptor-neurotransmitter dissociation followed by 
channel closing. We show that this lemma can be generalized and subsequently used for lifting 
most of the original constraints on the time constants. Thus we show that the Carnevale-Hines 
integration scheme for the integrate-and-fire model can be employed for simulating a much wider 
range of neuron and synapse type combinations than is apparent from the original treatment. 

PACS numbers: 87.19.11, 87.19.1g, 87.18. Sn 

Keywords: Event based, integrate-and-fire, neuron, synaptic current, Carnevale-Hines lemma, NEURON 



I. INTRODUCTION 

One of the most salient features of neurons is their ability to summate synaptic inputs arriving from other neurons 
and to respond with the generation of an action potential or spike when the membrane potential reaches a certain 
threshold value. After its generation, a spike will generally travel down the neurons axon to serve as an input to 
other cells, including muscles fibers and neurons. In its most basic form, spike generation is captured by the so-called 
integrate-and-fire model. This model was first conceived a hundred years ago by Lapicque [T] . Lapicque modeled the 
subthreshold behavior of the membrane potential as a capacitance in parallel with a resistor based on the electrical 
properties of the cell membrane. At that time, the spike generating mechanism was not known, and it was therefore 
only possible to give a phenomenological description of the process. On the basis of electrophysiological experiments, 
Lapicque assumed that when the membrane potential reached a threshold value, the cell would generate (fire) a spike 
and subsequently the membrane potential would be reset to resting level [U |2]. Integrate-and-fire models are still 
widely used today, both in simulations and for the analytical study of neural network dynamics. 

The integration scheme we analyze here was introduced by Carnevale and Hines [31 H] for the widely used NEURON 
simulation environment [5] and is event-based. In event-based models the synaptic coupling between neurons is 
mediated by events. Events are triggered by threshold crossings of the integrate-and-fire neurons and subsequently 
communicated to the postsynaptic cells. Postsynaptically these events initiate a change in the synapse, which in 
the most common case lead to excitation or inhibition of the postsynaptic cell. Support for event-based integration 
methods is available in several other scientific neural simulators ( e.g. NEST, XPP and Mvaspike) [B] which makes 
these simulators possible candidates for the implementation of the scheme discussed here. Besides the Carnevale-Hines 
scheme many other integration schemes are in use to simulate integrate-and-fire models. These schemes range from 
purely numerical integration schemes, such as Eulcr and Rungc-Kutta, to numerically exact calculations based on root- 
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finding algorithms for determining when the membrane potential crosses the spiking threshold [7J [HI E 113 HU H2I • The 
Carnevale-Hines scheme takes a middle ground between the two above mentioned extremes; it uses explicit knowledge 
of the exact solution to determine whether threshold crossing will occur, but avoids the expensive explicit calculation 
of the threshold passage time. Instead, the model employs the computationally cheaper Newton iteration to obtain a 
spike time estimate. For situations where presynaptic cells become active between two firing times and the order of 
firing of the cells is unknown (for example, due to mild external noise), we expect this scheme to be computational 
efficient. An analysis of computational efficiency, however, is outside the scope of this paper. 

The Carnevale-Hines scheme is developed to solve an integrate-and-fire model which includes excitatory synapse 
as exponentially decaying currents and inhibitory synapse as currents following a double exponential function. With 
these currents, it is in principle possible to describe the main class of excitatory synapses (characterized by AMPA- 
receptors) and the main class of inhibitory synapses (characterized by GAB A receptors). As mentioned before the 
Carnevale-Hines scheme uses Newton iteration to estimate the threshold crossing times, or phrased differently to find 
the events. In their proof of the correctness of the Newton iteration estimate Carnevale and Hincs used the following 
constraints on the time constants of the synapses and the membrane time constants: r^ecay, excitatory < Trise, inhibitory < 
Tdecay, inhibitory < ^membrane ED H] • These constraints imposed by the integration method lead to the loss of many 
physiological relevant realizations of the conceptual model. To elucidate the biological relevance of this point we will 
discuss the physiological parameter range in the next paragraph. 

In cortical areas of mammals, the excitatory AMPA currents have a fast rise time between 0.1 and 0.8 ms, followed 
by a fast decay of 1 to 3 ms (T3J H3] - In these areas, the inhibitory GABAergic currents have a rise time between 1 
and 2 ms |T5] and a decay time varying from about 5 ms to about 30 ms [T3]. Also the membrane time constants of 
different cortical neurons vary over a wide range, from close to 5 ms to well over over 40 ms [17j [T8j HH [20]. These 
physiological data show that the decay times for AMPA synapses are similar or larger than the GABA rise times and 
furthermore that the GABA decay times can be larger than the membrane time constant. However, as stated above 
the original treatment of the Carnevale-Hines integration scheme requires that the excitatory decay time should be 
smaller than the inhibitory rise time, and that the inhibitory decay time should be smaller than the membrane time 
constant. Consequently, the full physiological range of parameters as found in experiment is not accessible to the 
standard implementation of the model. These limitations became apparent to us during studies into the effects of 
GABA-receptor maturation on microcircuit processing, on which we will report elsewhere. This has lead us to the 
reexamination of the Carnevale-Hines integration scheme that we present here. 

The aim then of our analysis in this paper is to remove the unphysiological constraints on the time constants while 
keeping the strength of this event-based integration scheme. Our analysis will proceed in four stages. In the first stage 
we introduce the basics of the model. In the second stage we provide an analysis of the physics of the problem. In 
the third stage we offer our improved proof of the Carnevale-Hines lemma relaxing some of its original preconditions. 
In the fourth stage the actual analysis of the integration scheme is carried out using the generalized Carnevale-Hines 
lemma to abolish the unphysiological constraints on the decay times. 

II. MODEL AND INTEGRATION SCHEME 

In this section we will expand upon the short introduction of the event based integrate-and-fire model provided 
in the introduction. Subsequently we discuss the exact solution of the models subthreshold dynamics. Both in our 
treatment of the model and the exact solution we will in essence follow the original treatment [21 0], however our 
notation is slightly adjusted to account for the possibility of having synaptic subtypes. 

A. Event-based synaptic current driven integrate-and-fire model 

The model is event based, such models for neural networks are based on the assumption that communication between 
neurons is completely dependent on action potentials generated in the axo-somatic region and that as a result the 
time of occurence of the action potential contains all the information a single neuron communicates to the neurons it 
innervates. In the model used here an event is generated at the moment the neuron passes the firing threshold and 
axonal propagation and presynaptic delays arc accounted for by delivering it with an appropriate time delay typically 
in the order of a few milliseconds at a postsynaptic cell. The synapse in the postsynaptic cell responds then by 
generating a synaptic current. Physiologically, the time course of synaptic currents is determined by the association 
rate of the neurotransmitter to the receptor, the dissociation rate of the neurotransmitter from the receptor, the 
removal rate of neurotransmitter from the synaptic cleft and the driving synaptic reversal potential |21j . For current- 
based integrate-and-fire models, it is further assumed that the synaptic current follows the receptor opening, while 
the induced change in driving force due to a changing membrane potential is ignored. Here this approach is followed 
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and excitatory synapses and inhibitory synapses are included as membrane potential independent currents E v and 
J M , respectively. The subscripts v and /i in these expressions indicate the different excitatory and inhibitory synaptic 
subtypes. Because the currents follow the opening and closing of the receptors their time course is completely 
specified by the average time course of the open receptor-neurotransmittcr complex states. To describe these we 
use the variables e v and for the excitatory and inhibitory synapses, respectively. On the arrival of an event at 
an excitatory synapse the receptors in the excitatory synapse transfer to the open channel state instantaneously , 
reflecting a fast transient AMPA binding dynamics. In the model this is reflected by adding the weight of the synapse 
w e ,u > to e„. On the arrival of an event at an inhibitory synapse, however, we first get a fast increase of the amount 
of receptor-neurotransmitter complex in the closed state followed by a transition of these receptor-neurotransmitter 
complexes to the open state. This change in amount of receptor-neurotransmitter complex is modelled by adding the 
weight of the particular connection Wi „ < to an auxiliary variable j M describing the closed receptor-neurotransmitter 
complex state, the sign of the weight is used in this operation to indicate the inhibitory nature of the resulting current. 
After this event handling the variables e^j'^ and i M evolve according to the following differential equations linking the 
slowest relevant time scales of the kinetic models underlying receptor opening and closing of the AMPA and GABA 
receptor to a time development model for the synaptic currents: 

de v 1 

~jr = e u 

at t Cv 

a% _ _ J_ . 
It n 3,1 

din _J_- , 
dt ~ T ^ + a ^ 

(1) 

Next to the already introduced variables e u , and we also see the appearance of the time constants T e „ for the 
excitatory decay time, and Tj ,t; for the inhibitory rise and decay time. Next to these time constants which come 
with the dynamics we sketched, we see the appearance of a parameter where the appearance of tJ might have 
been anticipated, the parameter aj acts as normalization constant and is chosen such that an event induced change 
in jn by an amount uii , M results in a maximal change of w i fl in i„. For technical reasons and in line with the names 
chosen, we will assume that Tj < r» , which assumption basically restricts the mechanism underlying the double 
exponential current to mechanisms in which the slow time constant acts on % and not on j. Although this can be at 
odds with the actual biophysics, it poses no real constraint on this phenomenological model, which in its spike output 
is only sensitive to the shape of the inhibitory current which is insensitive to an interchange of Tj , Tj provided the 
normalization constants are adjusted accordingly . 

In the model, the membrane potential is represented by the variable m, and the resting potential and spike threshold 
are identified with m = and m = 1, respectively. When the membrane potential deviates from the resting potential, 
then in the absence of synaptic currents it decays back to resting potential and will do so with the membrane time 
constant r m . Furthermore, the excitatory I e — ^ v a ev e v and inhibitory li — X^ t a vV s y na Ptic currents act on the 
membrane potential. The constants a ev and ctj are normalization constants and are chosen such that an isolated 
instantaneous change in e v by an amount w e>v > induces a maximum depolarization of w eM and, similarly, such that 
an isolated instantaneous change in j v by an amount < induces a maximum hypcrpolarization of |H [S]. 
Alternatively, they can be chosen to normalize the charge transfer [S]. Putting everything together, the differential 
equation describing the membrane potential becomes: 

dm 1 x - x ^ . 

-rr = m + 2^ a ^ e " + Z^ a vV ( 2 ) 

(XL Tm. 



This differential equation together with the accompanying differential equations for the synapses has an exact 
solution for the subthreshold behavior. These differential equations are linear in the currents and therefore the 
solution is completely analogous with the case for a single excitatory and a single inhibitory current [4, 5 . Using the 
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reciprocals of the time constants k x = 1/t x we can write the exact solution as follows: 



e„(t) = e,,oe- fe - (t -* o) , (3) 

3,(t) = jpfie-W-^, (4) 

- V,oe- fe ^ (t -* o) , 

+J M A>~ Mt ~ to) -e~ fcj '' (t ~ to) ) (5) 

m(t) - m e~ fcm( *~* o) 

+ S e ".o 6 ^( e "* m( *" to) - e" fee " ( *"* o) ) 

+ Ev.o^(e- fcm(t -* o) -e- fc ' (t -* o) ) 

_ fey - k ra ^ c - km (t-t ) _ g-fe^ (t-t ) )]_ 



(6) 



In this expression to refers to the last time preceding f at which m, e^, j M and i M were evaluated, the values of those 
variables at to are denoted by mo, ev,o,jnfii V,o- The constants a e „, and a v used in the differential equations are 
absorbed into new constants 6 together with part of the k x , i.e. b e — a e /(k e — k m ),bi — di/(ki — k m ),bj — a,j/(kj — fcj). 
The challenge entailed in this subthreshold behavior is to extract from it the point where the membrane potential 
crosses threshold, which is the topic of the next section. 



B. Carnevale-Hines integration scheme 



The Carnevale-Hines integration scheme is developed to solve the problem of finding threshold passage times from 
the exact solution of the model presented in the previous subsection. The scheme cycles through the following steps: 
in the first step an event arrives at the neuron, in the second step the actual membrane potential and synaptic 
currents at the arrival time of the event are calculated from the exact solution, in the third step a check takes place on 
threshold crossing and in the fourth step the actual event is handled by updating synaptic currents and calculating a 
new threshold crossing estimate based on these currents, in the fifth step a self-event which will arrive at the estimated 
threshold crossing time is generated. Let us examine steps three and four in slightly more detail. In the third step a 
test takes place to establish whether the membrane potential m is close to the threshold 9—1 (i.e. m > 9 — e, with e 
an arbitrary number satisfying < e << 9). If sufficiently close to threshold then it is assumed that the membrane 
potential reached threshold and the membrane potential m is reset to zero and an event is sent to the synapse of the 
cells innervated by this cell. In the fourth step the event is handled. If the event was a synaptic activation then the 
weight of the synapse is, inline with the model description above, added to either an e„ or an j M . If the event was a 
self-event related to an estimated threshold crossing time then the currents obtained from the exact solution are kept 
unaltered. After this update of the synaptic currents an estimate is made about when to evaluate the subthreshold 
solution again. This estimate is based on the time derivative of the membrane potential m! which is equal to the total 
of the synaptic and leak currents. If the current is not depolarizing ml < 0, then the evaluation of the membrane 
potential is postponed until the arrival of a new synaptic event. If the current is depolarizing m! > 0, a future spike 
time is estimated t e using Newton iteration t e = (1 — m) /m' + 1 . This estimated time is used in step five to generate 
another type of event, a self-event that is sent by the neuron to itself and indicates the latest time to which evaluation 
of the membrane potential for threshold detection can be postponed. For this approach to work it is necessary that 
the estimated threshold crossing time is before the actual threshold crossing, so that the exact solution can be used 
to evaluate the membrane potential at that time and we can simply start our event handling again to detect whether 
we actually crossed threshold. 

The essential assumption of this mechanism is that every Newton iteration step is underestimating the threshold 
passage time, which allows it to approach the threshold crossing with several iterations without the risk of overesti- 
mating it. The discussion in the next section and the subsequent mathematical proof focus on showing that provided 
the excitatory synaptic currents decay faster than the inhibitory currents (r ev < t v for all combinations of v and \x) 
the Newton iteration step is underestimating threshold passage time. 
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III. THRESHOLD PASSAGE TIME IS UNDERESTIMATED BY NEWTON ITERATION 

The formal proof we present in this section is largely analogous to the proof in the original analysis [3l|5], except 
that we base it on a generalization of the underlying Carnevale-Hines lemma, for which we give a proof here. This 
generalization allows us to lift the unphysiological constraints in the original analysis. The actual proof consist of 
two parts: the first part shows that when at to the derivative to' = dm/dt < 0, the membrane potential will stay 
below threshold at least until a new event arrives; the second part shows that when at to the derivative m' satisfies 
mf > 0, the Newton iteration formula t e = (l — m)/m' + t underestimates threshold passage time. Before presenting 
the formal proof based on the exact solution for m(t), we analyze the physics, which provides better insight into the 
actual underlying mechanisms. 



To answer the question when Newton iteration underestimates threshold passage time we need to look at the different 
components contributing to the membrane potential derivative to' — dm/dt. We have three kinds of currents: a leak 
current, which acts in the direction of the resting membrane potential; excitatory synaptic currents, which depolarize 
the membrane; and inhibitory currents, which hyperpolarize the membrane. If the Newton iteration estimate of 
threshold passage time t e 



is required to underestimate threshold passage time, then that puts limitations on the possible time course of the 
currents contributing to to'. If the membrane potential does not cross threshold 6 between to and t e , then it satisfies 
the inequality 



Using simple arguments, we can derive conditions on the time constants for which the sum of synaptic and leak 
currents is decreasing. From the inequality above we can immediately see that if during t < t < t e the total current 
shows no growth, i.e. m'(t) < m' , then no threshold passage will take place in this time interval. Therefore, under the 
stronger assumption that the total current is decreasing the Newton iteration will underestimate threshold passage 
time. 

Let us first analyse the situation in which no inhibitory currents are present (for example the top black lines in 
figure flj. We know that the excitatory currents are decaying, i.e. their contribution to to' will reduce over time, and 
the only uncertainty is in the leak currents. If to' < the excitatory synaptic currents cannot overcome the leak 
current at the present membrane potential, and therefore after decaying further they are definitely unable to do so. 
On the other hand, while m' > the leak current is growing and therefore reduces mf, and we see that m! < m 
until threshold is reached or until ml reverses sign and the membrane potential moves away from threshold never to 
reach it. From this we see that either the threshold will never be reached and every finite estimate is underestimating 
threshold passage time or to' satisfies the inequality m'(t) < m' upto reaching the threshold and we know then that 
no threshold passage takes place before the estimated time t e . 

Now, let us examine the situation in the presence of double exponential inhibitory currents; again examples are 
given in figure [I] where grey lines indicate cases where the inhibitory currents are growing and black lines cases 
where there are only decaying inhibitory currents. The first observation is that a growth of the inhibitory current 
leads to reduction of to', so adding the growth component associated with j M will only strengthen the arguments 
we can obtain after assuming it to be equal to zero. So assuming j M = and m' Q > and mo > we know that 
the individual synaptic currents are decaying towards 0. If, however, the excitatory contribution would decay slowly 
while the inhibition would decay fast, then the actual resulting synaptic current might be a growing depolarizing 
current and we might overestimate firing time. It is this consideration that leads to our only real constriction on the 
time constants: the decay times for excitatory synaptic currents should be faster than those for inhibitory synaptic 
currents. If, on the other hand, the inhibitory currents decay slower than the excitatory currents, we know that the 
total synaptic current decreases faster than expected on the basis of the excitatory decay time constants alone and 
might even become hyperpolarizing. Again we can see that while to' > the leak current is growing and therefore 
also reduces to'. Taken together, we find that to' < m' until threshold is reached or until m' reverses sign and the 
membrane potential moves away from threshold never to reach it. So under the assumption that inhibitory synaptic 
currents decay slower than excitatory synaptic currents, the condition to' < m' is fulfilled and we know that the 
estimated threshold passage time t e based on Newton iteration underestimates threshold passage time. If contrary to 



A. When will Newton iteration underestimate threshold passage time? A physical analysis 



t e — {9 — mo) / to + to f° r m 'o > 



(7) 




(8) 
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FIG. 1: Illustration of the physical analysis. The figures show the derivative of the membrane potential, m', against time. 
In a., the whole time course is shown, and in b. only the initial part. If mo,m' > then m! will stay below m' upto the 
Newton estimated threshold passage time indicated by the vertical red bar. The top black line represents the case were initially 
all synaptic currents are excitatory and the inhibitory synaptic current is zero for all times; it is shown upto the point were 
threshold is actually reached. The lower lying black lines are at larger initial inhibitory (and hence excitatory) synaptic currents 
(m'i = a%i — —0.20, —0.15, —0.10, —0.05) but without the growth component (j — 0). The accompanying grey lines show the 
effect of adding the growth component (j = —0.025, —0.050, —0.075, —0.100). Parameters used: r m = 30, r e = 3, Tj = 9, Tj — 2, 
m = 0.5, m'o = 0.25. 

our earlier assumption too < 0, then we have two separate cases: the first case where the excitatory synaptic current 
is larger than the inhibitory synaptic current, in which situation all the arguments above apply; and a second case in 
which we have a dynamics dominated by the leak current. For this leak-dominated phase the argument above does 
not apply, but because the leak current acts towards the resting potential and not towards the threshold, it is clear 
that no threshold passage will take place and any estimate will underestimate threshold passage time. 

From these arguments we also see why it is difficult to include NMD A like currents, which are best modeled 
by an excitatory double exponential current. Their growing excitatory contribution accelerates the rate at which 
the membrane potential approaches the threshold after activation of the synapse, causing the Newton iteration to 
overestimate the threshold passage time. 

B. Generalized Carnevale-Hines Lemma 

The exact expression of the membrane potential is built from a large number of double exponential functions. 
Establishing upper (lower) bounds on such a sum of double exponentials is strongly simplified by the Carnevale-Hines 
lemma, which allow us to replace one double exponential with another that for all times is larger (smaller). The 
lemma follows from the following corollary, which gives us the monotonic development of double exponentials when 
viewed as a function of one of the decay times: 

Corollary III.l. For i,/i,Ae$ the functions 



/a(m,*) 



e -Ai _ g-^t 
)JL — A 



(9) 



are defined for /i ^ A and by including the limits lim^x f\(fi,t) into f\ it can be extended to a function continuous 
in \x. For fixed t the extended function f\ is a monotonically decreasing continuous function of \x. 

Proof. We start with showing that the function f\ can be extended to a continuous function of n by adding the point 
fi = A. The numerator and denominator used in the definition of f\ are at /i = A, and their derivatives with respect 
to fi exists; furthermore, the derivative of the denominator is equal to 1 and therefore non-zero. From this we can 
see that the preconditions of l'Hopital's rule are satisfied. From l'Hopital's rule we know that the limit exist and we 
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find lim M ^A f\(/i,t) — te~ M *, i.e. the alpha function. After extending f\ with this limit, continuity in fi follows from 
observing that l'Hopital's rule is based on the fact that the left and right limits are equal and hence no discontinuity 
occurs at fj, — A in the extended f\. 

The next step is to prove that the derivative of the function fx with respect to /x is negative a.e. for t ^ 0. This 
derivative is given by: 



d[i (n — xy 



- A) + 1 - e^*) . (10) 



On the right hand side of this equation the first factor is clearly positive, and therefore we need to prove that the 
second factor is negative a.e. to show that the function is monotonically decreasing. To show that the second factor 
is negative a.e., we prove that it has a non-positive maximum at A = fi. If we examine the derivative of the second 
factor 

d(fc-A)t + l-e("-A)') _ _ t 



dfi 



= t-te (t *- X)t (11) 



we find that it is for /i = A, positive for \x < X and negative for /i > A, showing that there is indeed a maximum 
(t(fj. — A) + 1 — e' M_A - ) *) = in the second factor at A = /i. As a result we find that df\(fi, t) /dfx < for fi ^ X and the 
only step left is to prove that df\(fi,t)/d[i is continuous. 

The expression for df\(/i,t)/d^i is indeterminate at \x — A, but the preconditions for l'Hopital's rule are satisfied, 
indicating that this derivative is continuous at this point. The values of the df\(fj,,t)/dii at /i — X are most easily 
calculated by inserting the Taylor expansion for t(fi — A) + 1 — e^ - ^ 1 into equation 
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^M = _ e -/, t g (/x _ A) »-2*" (12) 

d ^ n=2 U - 

The indeterminacy at fi = X is now canceled and we find 

_ (13) 

This expression is negative and thus the derivative df\(fj,, t)/dfi is negative everywhere. □ 
Lemma III. 2 (Generalized Carnevale-Hines Lemma). If fi 2 > fJ-i and MiiM2 7^ A then 

g— At _ e — mt e~ A * — e~^ 2 ' 



Mi - A fi 2 — X 

Proof. Equality follows from choosing t = for which we have /a(m, 0) = for all values of /i. For t ^ the inequality 
is a direct consequence of the corollary. □ 



C. Movement away from threshold implies threshold will never be reached 

We start this analysis from the exact solution of the membrane potential m(t) given in equation [6] The purpose is 
to show that the exact solution has an upperbound given by the line m e (t) = m + rri' (t — to). Instead of showing 
that this is true because m'(t) < m' as in the physical analysis, we will now show it on the basis of the exact solution 
itself. 

As before, the terms describing growth of the inhibition, i.e. those terms related to j^, can be discarded, but now 
we use the Carnevale-Hines lemma to achieve this. From the lemma we find that for r,- < t, the factors, 

e -k m (t-t ) _ e -h ilM {t-t ) _ k ^ Z km L-k m (t-t ) _ e -k jlt {t-t )\ > q ^ 

found in the exact solution for m(t) (equationj6j) are all positive. The factors multiplying these expression, b ifi and 
bj , are also positive, but the factors j Mj o are negative. The terms, therefore, in which these appear are negative and 
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we obtain an upperbound for m by dropping these terms: 

m(t) < m a e~ km{t ~ to) 



J2 { ^K ( e fe " 



(*-*o) _ e -fci M (t-*o) 



(16) 



In the next step we use the assumption that t Cv < for all combinations of [i and v to be able to apply the 
Carnevale-Hincs lemma. Using an arbitrary excitatory decay time r ey we replace the dependence on Tj in the i^o 
related terms with a r e , dependence. If we fix r e , by choosing the largest excitatory decay time for it, we can use 
the lemma a second time and replace the r ev dependence in the e v $ related terms with a r e , dependence as well: 

m(t) < m e- fcm( *- to) 



' \ e -fc m (t-«o) _ e -*v(*-*o) 



Now we can use m' — —k m m + ^ v a &v e v + ai^i^ to obtain: 

m(t) < m e^ km{t ~ ta) 

+m h \ (e- k ^-^ - e- k <»> {t - to) ) 

+m' 1 ( e -*»(«-*o) _ g-*.,, (t-to)j ( 18 ) 

We assumed that m' < 0, and because the other factors multiplying it are all positive, we can remove the associated 
term from the inequality, so that we obtain, 

m{t) < m e- fcm(t - to) 

+m k ™ ( e -fcm(«-*0 - e -*.^(*-*o)) (19) 

This expression is now fully equivalent to the one found in the original treatment [H [5] , where it is shown that the 
derivative with respect to time of the term multiplying too is negative; and because it is 1 at t = to , it will be smaller 
than 1 at larger times while it is also positive. Because Too < 1 this shows that m(t) stays between — oo and 1 and 
hence below threshold. 

D. Movement toward threshold shows slowdown 



If to' > the step from equation 18 to equation 19 is not allowed, but using our corollary we can replace the last 



term in equation 18 with the alpha function belonging to the slowest decay time out of r m and T e y: 

m(t) < TO e"' £m(t " to) 

+m k ™ ( e -M«-*0 _ e-^C-'o)) 

+mf(t - to)e- mm{k ™< k ^> )(*-*<>) (20) 

By the same argument as used before, the first two terms combined are smaller than or equal to too; also 
e -m%n(k m ,ke v ,)(t-t ) - 1S sma ]j ei . than 1, so taking this together with too, (t — to) > we obtain: 

m(t) < m + m'(t-t ) (21) 
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When mo < 0, the second term in the right-hand side of equation 20 is negative and we can drop it from the inequality, 
which is the correct thing to do if T e y < r m . If r e y > r m we can reorganize the first two terms in such away that 
we can replace it with e~ fcc >"'( t ~ t °) : 

m(t) < m (e-^-V + km ( e - fc m(t-«o) _ e -^(t-*o)A 

+m'(t - t )e- mmikm - k "^ ){t - to) 

( ke r e -k m (t-t ) _ e -*v (*-*oh + e -K : „,(t-t )\ 

V ke' v k m J 



m 



+m'(t - t ) e - min ( fc "" fc «,' )('-'") 



(22) 



From which we see that we have a free choice on whether we want to keep the factor e k "- 4 °) or e fcm ^* to ' from 
the first two terms. It will be convenient to make a choice which fits the last term 

m(t) < (m + m '(t-to))e~ mmikm ^' ){t ' to) (23) 

This expression shows that in the area where the linear extrapolation is below zero we have no guarantee that m(t) 
is under the linear extrapolation, but the actual upper bound is below zero and no threshold passage occurs. When 
the linear extrapolation crossed zero, we again find m(f) < mo + m'(t — to), and we see that Newton iteration is 
underestimating the threshold crossing time. This seemingly strange behavior of the upperbound is not indicating 
a magical zero crossing of m(t) at the point where m + m'(t — to) changes sign, but reflects our lack of knowledge 
about the system. If the system was dominated by excitatory currents during the whole period, then we know from 
our physical analysis that the linear extrapolation was above the real curve all the time. If however the system was 
dominated by an unmasked leak current, then the real curve went above the linear extrapolation, but because the 
whole movement is towards resting potential and will never cross the resting potential, the linear extrapolation has 
to cross the membrane potential before passing through zero. 



E. Summary 



Let us summarize the behavior of the membrane potential in this model as we analyzed it both by the integrated- 
current based argument and by the Carnavale-Hines lemma. When excitatory synaptic currents dominate (see equa- 
tions 21|23 1 , as in the early stages of the examples given in figure [2] the membrane potential increases but always 
stays below the linear extrapolation . When inhibitory synaptic currents dominate (see equation 19 1, the membrane 



potential moves away from threshold, and if they are sufficiently strong, as in the lower curves shown in figure [2] 
they lead to reversal of the polarity of the membrane. When in the last part of these curves leak currents dominate 
(see equations 19|23 1, the membrane potential moves toward resting potential. If the latter takes place at negative 
membrane potentials, the membrane potential can exceed linear extrapolation, because decay of inhibitory synaptic 
current can quickly unmask an upward leak current. Although in this situation the linear extrapolation is exceeded, 
no overestimation of threshold passage time will occur because no threshold passage takes place in the leak dominated 
phase. In conclusion, we established that either no threshold crossing will take place or that threshold crossing takes 
place after the Newton estimate for threshold crossing time. 



IV. DISCUSSION 



We have shown that the non-physiological constraints on the time constants in the event-based integrate-and-fire 
Carncvale- nines integration scheme can be relaxed. The Newton iteration step is underestimating threshold passage 
time provided the slowest decaying excitatory synaptic current (i/ a ) decays faster than the fastest decaying inhibitory 
current ((J>f), i.e. we have (T e ^ s < ) . This makes the model applicable to a much wider range of neuron and 
synapse combinations in the nervous system. Although we here analyzed a current-based model, the argument about 
threshold passage time, which is based on our physical analysis, can with slight modifications also be applied to 
conductance-based integrate-and-fire models. Conductance-based models have synaptic reversal potentials, but these 
do not weaken our argument, as can be seen as follows. On approaching the reversal potential, the driving force for 
the excitatory current is reduced. The excitatory currents will therefore be reduced further than would be expected 
from channel kinetics alone. The driving force of the inhibitory currents is enhanced on approaching the threshold 
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potential. The inhibitory currents will therefore be less reduced than would be expected from channel kinetics alone. 
We see from this that including reversal potentials will only strengthen the threshold passage time argument. The 
only ingredient missing, then, to apply the Carncvale-Hines integration scheme to conductance-based integrate-and- 
fire models is the lack of an exact solution for the membrane potential or another numerically efficient way to calculate 
the membrane potential after a time step of arbitrary size. Within the Carncvalc-Hincs integration scheme based on 
Newton iteration slowly growing excitatory currents (like NMDA currents) cannot be included. We can however, 
sometimes, replace the Newton iteration by an higher order extrapolation scheme. Lets for example take a synaptic 
model for NMDA currents using the same differential equations as are used for the inhibitory current (admittedly 
ignoring the magnesium block, which is not within reach of our treatment here). We now take i,j to represent the 
double exponential excitatory currents and auxiliary variable, respectively, and leave out the inhibitory currents. To 
do this we only need to change the sign of j, so we assume j > 0. If we now further assume that T{ > r m then the 
exact solution has an upperbound. 

m(t) <m + m't + j ^ - (t - t ) 2 (24) 

This quadratic upperbound for the exact solution was found using appropriate a-functions as upperbounds for the 
double exponentials, i.e. by repetitively using the generalized Carnevale-Hines lemma. We again obtain a threshold 
crossing estimate if we determine where this upperbound crosses threshold. This new estimate does not solve the 
problem of unmasking yet, but we expect that unmasking can be dealt with in a similar way leading to extra terms 
quadratic terms in the upperboundary. The reason is that only already existing decaying currents can be unmasked, 
those currents lead to a double exponential contribution to the membrane potential development and putting an 
upperbound on a combination of two double exponentials involves two times in succession the replacement of a double 
exponential by an a-function leading to a (f — to) 2 term. If this reasoning is correct then in even more general 
physiological circumstances, i.e. T decaVtAMPA < T decay , G ABA < T decay ^ NMDA a threshold crossing estimate can be 
found by taking the positive root of a quadratic polynomial. 




FIG. 2: Illustration of the Carnevale-Hines lemma based argument. The figures show the membrane potential, m, against 
time. In a, the whole time course is shown, and in b only the initial part. If mo, m' > then m will stay below mo + m' t up 
to the estimated threshold passage time indicated by the point were the tangent at to (red line) crosses threshold 9 (dashed 
red line). The top black line represents the case were initially all synaptic currents are excitatory and the inhibitory synaptic 
current is zero for all times; it is shown up to the point were threshold is reached. The lower lying black lines are at larger 
initial inhibitory synaptic currents (mj = 041 — —0.20, —0.15, —0.10, —0.05) but without the growth component (j — 0). The 
accompanying grey lines show the effect of adding the growth component (j = —0.025, —0.050, —0.075, —0.100). Parameters 
used: r m = 30, r e = 3, n = 9, Tj = 2, m = 0.5, m' = 0.25. 
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